探究某个基因的外显子覆盖度情况【直播】我的基因组87

一般情况下，我们得到了测序reads在基因组的比对情况文件bam格式的，里面的信息非常多，如果我想特定的查看某个基因的情况，那么我们可以选择IGV等可视化工具，但它并不是万能的，因为即使是一个基因，它也会有多个转录本，多个外显子。以前我写过批量IGV截图(点击直达)，但是大部分基因的长度都超过了37Kb，超过了IGV的窗口浏览限制。而且我们也不需要知道该基因上面比对成功的所有reads信息，太复杂了，我们只需要知道基因上面各个部位的测序深度即可，而且基因上面比较重要的就是外显子了，被一个个内含子隔离开来。

所以，我们可以画它的外显子覆盖图，下面是一个例子：横坐标是外显子的长度，纵坐标是测序深度，每一个小图都是一个外显子

DMD外显子覆盖度情况

根据这个图，我们就可以很明显的看出，DMD基因NM_000109转录本的1，10-17号外显子缺失，用IGV一个个的看这些外显子区域，是同样的结果！可能是芯片捕获不到，也可能是样本本身变异，造成的大片段缺失。但是这个图的信息就非常有用！(当然，这个肯定不是我的啦，我很正常的哦)

PS：请忽略上图的外显子不是按照数值的大小排序，只是因为当初我对ggplot还不是很熟悉，不知道调整factor就可以改变出图的顺序。

那么，我们该如何画这样的图呢？

首先，我们需要找到需要探究的基因的全部转录信息，及外显子信息！

这里的hg19_refGene.txt我直接从annovar的数据目录拿到的。可以看到一个基因有多个转录本，每个转录本的外显子个数不一样。如下：

那么，我们根据这个信息，就可以判断该基因的每个外显子的起始终止位点啦！

然后用samtools的depth命令去找这个基因的全部片段的测序深度信息

最后再格式化成下面的三列数据

1. 1    226 exon:43

2. 2    235 exon:43

3. 3    246 exon:43

4. 4    254 exon:43

5. 5    258 exon:43

6. 6    262 exon:43

7. 7    277 exon:43

8. 8    286 exon:43

9. 9    298 exon:43

10. 10    319 exon:43

• 第一列是该外显子的起始终止坐标，从1到该外显子的长度。每切换一个外显子，坐标从1开始记录。一般的外显子长度在200bp左右，也有一些超长的外显子5kb及以上长度。

• 第二列是该外显子在该坐标的测序深度，通过samtools的depth命令得到

• 最后一列是该外显子的标记，从exon:79一直倒推到exon:1，因为该基因在染色体的负链，所以外显子顺序是反着的！这一列信息用来把外显子在ggplot上面分面！

最后根据这个txt文档，用R语言，很容易就画出上面那样的图片了！R语言程序是：

1. if("ggplot2" %in% rownames(installed.packages()) == FALSE) {install.packages("ggplot2")}

2. library(ggplot2)

3. args <- commandArgs(trailingOnly = TRUE)

4. file=args[1]

5. outpng=sub(".txt",".png",file)

6. dat=read.table(file)

7. names(dat)=c('pos','depth','exon')

8. new_order=paste0('exon:',1:length(unique(dat[,3])))

9. dat[,3] <- factor(dat[,3] ,new_order )

10. png(outpng,res=120,width = 1080, height = 1080)

11. p=ggplot(data=dat,aes(x=pos,y=depth,color=exon))+geom_line()

12. p=p+facet_wrap(~exon,scales="free_x")

13. p=p+theme(legend.position='none')

14. print(p)

15. dev.off()

这个是修正之后的R代码，外显子的顺序ok了。

比如下面这个捕获测序的，可以很明显的看到外显子区域测序深度超级高！但是第一个和最后一个外显子很明显测序深度不足，有可能是设计这个panel的人认为第一个和最后一个外显子是UTR区域，不予捕获。

用这个代码画一下我的WGS数据

1. bam=$1

2. file=$(basename $bam )

3. sample=${file%%.*}

4. echo $sample

5. mkdir -p  exon_png_$sample

6. genes=(DMD TP53 BRCA1 BRCA2 ALK ROS1 EGFR MET BRAF FGFR1 RET KRAS)

7. for gene in ${genes[@]}

8. do

9. perl ~/biosoft/exoncoverage/bin/cal.pl  -r ~/biosoft/exoncoverage/db/hg19_refGene.txt \

10. -b $bam  -S $gene

11. done

12. mv *.png  exon_png_$sample

13. rm *.txt

可以看到，我只是挑选了部分基因来画图，大部分都挺好的， 全部覆盖到了，但是有一个BRCA1的转录本的其中一个外显子看起来是缺失了。如下：

可以看到我的WGS测序深度在50~60之间，这很正常！因为我就花了这些钱而已。

不过BRCA1的第2，8外显子几乎没有reads覆盖，这个很诡异，但是因为这两个外显子长度小于100bp，很有可能是NGS测序技术的限制，所以我倒是没那么担心。

如果大家有人类的WGS测序数据的，可以用我的方法做一个类似的图，跟我交流一下，我的邮箱是 [e-mail]  希望你可以给我一点信心，我很正常。

猜你喜欢

菜鸟入门

数据分析

编程实践

直播基因组分析

生信技能树

点击下面的阅读原文可以拿到本次分析的所有代码！